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Abstract. In this paper, we study a time discrete scheme for the initial value problem of 
the ES-BGK kinetic equation. Numerically solving these equations are challenging due to 
the nonlinear stiff collision (source) terms induced by small mean free or relaxation time. 
We study an implicit-explicit (IMEX) time discretization in which the convection is explicit 
while the relaxation term is implicit to overcome the stiffness. We first show how the 
implicit relaxation can be solved explicitly, and then prove asymptotically that this time 
discretization drives the density distribution toward the local Maxwellian when the mean 
free time goes to zero while the numerical time step is held fixed. This naturally imposes 
an asymptotic-preserving scheme in the Euler limit. The scheme so designed does not need 
any nonlinear iterative solver for the implicit relaxation term. Moreover, it can capture the 
macroscopic fluid dynamic (Euler) limit even if the small scale determined by the Knudsen 
number is not numerically resolved. We also show that it is consistent to the compressible 
Navier-Stokes equations if the viscosity and heat conductivity are numerically resolved. 
Several numerical examples, in both one and two space dimensions, are used to demonstrate 
the desired behavior of this scheme. 
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1. Introduction 

When gas is in thermal non-equilibrium, which are encountered frequently in hypersonic 
flows, vehicles at high altitudes and flows expanding into vacuum, the macroscopic constitu- 
tive laws based on the continuum hypothesis tend to breakdown. A critical parameter that 
characterizes the rarifiedness of the gas is the Knudsen number (e = A/L), where A is the 
average distance traveled by the molecules between collisions, or the mean free path, and 
L is the characteristic length scale. When the flow gradients are large, such as in shock or 
boundary layers, continuum fluid dynamics equations are not adequate, and one needs to 
use a kinetic equation. The fundamental kinetic equation for rarefied gas is the Boltzmann 
equation, 

(1.1) ! + „. Vl / = i £>(/,, 

which governs the evolution of the density f(t,x,v) of monoatomic particles in the phase, 
where (x,v) £ M. dx x M. dv . Boltzmann's collision operator has the fundamental properties of 
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conserving mass, momentum and energy: at the formal level 



Q(f,f)<f>(v)dv = 0, <j>(v) = l,v,\v\ 



Moreover, the equilibrium is the local Maxwellian distribution 



(1.2) 



M[f](v) 



P 



(2irT)*»/ 2 



exp 



\u — v\ 
2T 



where p, u, T are the density, macroscopic velocity and temperature of the gas, defined by 



(1.3)p 
and 



f(v)dv 



M[f](v), u = - 
P 



v f(v) dv = - 
P 



vM[f](v) dv, 



(1.4) r =T~ / \u-v\ 2 f{v)dv = ^— f \u-v\ 2 M[f]{v)dv. 

dvP jR d v d v p J^d v 

The Boltzmann equation is closely related to the Navier-Stokes system which governs the 
evolution of macroscopic density, momentum and energy in the continuum regime: 

dp 



(1.5) 



at 



+ div x (pu) = 0, 



dpu 

~~dT 



+ diva (pu®u + pi) = ediv x [p a(u)], 



dE 

— + div x ((£ + p)u) 



ediv x {po~{u)u + kV x T). 
where p is the pressure, E represents the total energy 



E 



1 



d, 



z i — u m 

2 pu + Y pT > 



and I is the identity matrix. Moreover, the tensor cr(u) denotes the strain rate tensor given 
by 

2 

a(u) = (V x u+ (Vxuf) - —div x ul. 

cl v 

These equations constitute a system of 2 + d v equations in 3 + d v unknowns. The pressure is 
related to the internal energy by the constitutive relation for a polytropic gas 

1 



P 



( 7 - 1) IE 



■ pu 



where the polytropic constant 7 = (d v + 2)/d v represents the ratio between specific heat at 
constant pressure and at constant volume, thus yielding p = pT while the viscosity p = p(T) 
and the thermal conductivity k = k(T) are defined according to the linearized Boltzmann 
operator with respect to the local Maxwellian [2]. 

Since the quadratic collision operator Q(f) has a rather complex form, simpler models have 
been introduced. The main requirement is to build models that have the right conservations, 
entropy condition described by the //-theorem, and have the fluid dynamics (Euler and 
Navier-Stokes) limits with the correct transport coefficients. The simplest model is the so- 
called BGK model introduced by Bhatnagar, Gross and Krook [4j. It is based on relaxation 
towards the local Maxwellian 

(1.6) Q(f) = - e {M[f] - /), 

where r depends on macroscopic quantities p, u and T. 

This model conserves mass, momentum and total energy, and has the correct Euler limit 
when e — > 0. But in the Chapman-Enskog expansion, the transport coefficients, that is p and 
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K obtained at the Navier-Stokes level are not satisfactory. In particular, the Prandtl number 
defined by 

7 — 1 K 

which relates the viscosity to the heat conductivity, is equal to 1, whereas for most gases, we 
have Pr < 1. For instance, the hard-sphere model for a monoatomic gas (7 = 5/3) in the 
Boltzmann equation leads to a Prandtl number very close to 2/3. 

One model, proposed by Holway [19J, has the desired property of having the correct con- 
servation laws, yields the Navier-Stokes approximation via the Chapman-Enskog expansion 
with a Prandtl number less than one, and yet is endowed with the entropy condition pQ. See 
also 0|9]. This model is referred to as the ellipsoidal statistical model (ES-BGK), where the 
Maxwellian Ai [/] in the relaxation term of (|1.6p is replaced by an anisotropic Gaussian Q [/] . 
In order to introduce the Gaussian model, we need further notations. Define the opposite of 
the stress tensor 

(1.7) Q(t,x) = - [ (v-u)®(v-u)f(t,x,v)dv. 

P jR d v 

Therefore the translational temperature is related to the T = tr(Q)/d v . We finally introduce 
the corrected tensor 

(1.8) T(t,x) = [(l-v)TI + i/0] (t,x), 

which can be viewed as a linear combination of the initial stress tensor and of the isotropic 
stress tensor TI developed by a Maxwellian distribution. 

The Gaussian model introduces a corrected BGK collision operator by replacing the local 
equilibrium Maxwellian by the Gaussian G[f] defined by 

(v — u) T^ 1 (v — l 



Q[f] 



p 



cxp 



■ N /det(27r70 

Thus, the corresponding collision operator is now 
(1-9) Q(f) = - e {Q[f] 



where r depends on p, u and T, the parameter 
the Prandtl number through the formula [lj 



T/2<i/<lis used to modify the value of 



< Pr 



1 



3 ~ l-v 
It first follows from the above definitions that 



< +00 for v G [-1/2, 1]. 



f(v) dv 



v f(v) dv 



g[f](v)dv = P , 



vQ[f}{v) dv = pu, 



f(v) dv 



Q[f]{v)dv = E 



and 



v 2 JM. dv ™ 

(v — u) <8) (v — u) f{v) dv = pQ, 

(v — u) ® (y — u) Q[f] dv = pT. 



This implies that this collision operator does indeed conserve mass, momentum and energy 
as imposed. The collision frequency, v, involves the Prandtl number Pr as a free parameter. 
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This allows the ES-BGK collision model to reproduce transport coefficients, viscosity and 
thermal conductivity, in the Chapman-Enslog expansion [2j [5] , recovering the Navier-Stokes 
equations density p, momentum pu and temperature T, with the correct Prandtl number. 

In this paper we study a temporally implicit-explicit (IMEX) discretization of the ES- 
BGK model. The advantage of such a time discretization is that it is uniformly stable with 
respect to the small Knudsen number, thus removing the stiffness of the relaxation term, yet 
the implicit relaxation term can be solved explicitly, thanks to the special structure of the 
relaxation term. Although such a property was realized for the classical BGK operator |12j . 
the ES-BGK operator is different, and we realized that one has to compute to the higher 
moment in order to evaluate the implicit Gaussian distribution explicitly. Furthermore, we 
show that this time discretization is asymptotic-preserving [20], an important property for 
the scheme to be robust in the fluid dynamic regime, allowing it to capture the fluid dynamic 
behavior without resolving numerically the small Knudsen number. We further show that 
this discretization is consistent to the compressible Navier-Stokes equations with the correct 
Prandtl number, if the viscosity and heat flux terms are suitably resolved numerically. We 
then validate this numerical method by presenting numerical results based on this scheme, 
and compare them with the solutions of the Boltzmann equation and the corresponding 
Navier-Stokes equations, in both one and two space dimensions. 

2. An Asymptotic Preserving scheme to the ES-BGK equation 

Past progress on developing robust numerical schemes for kinetic equations that also 
work in the fluid regimes has been guided by the fluid dynamic limit, in the framework 
of asymptotic-preserving (AP) scheme. As summarized by Jin |20j . a scheme for the kinetic 
equation is AP if 

• it preserves the discrete analogy of the Chapman-Enskog expansion, namely, it is a 
suitable scheme for the kinetic equation, yet, when holding the mesh size and time 
step fixed and letting the Knudsen number go to zero, the scheme becomes a suitable 
scheme for the limiting Euler equations 

• implicit collision terms can be implemented explicitly, or at least more efficiently than 
using the Newton type solvers for nonlinear algebraic systems. 

We now introduce the time discretization for the ES-BGK equation (jl.ip . (jl.9p 



(2.1) 



+ v V x f = ~ (£[/] - /), xenc R dx , v e R d \ 



k f(0,x,v) = f (x,v), xen,v£R d \ 

where r depends on p, u and T. 

The time discretization is an implicit-explicit (IMEX) scheme. Since the convection term 
in (|2.ip is not stiff, we will treat it explicitly. The source terms on the right hand side of 
(|2.ip will be handled using an implicit solver. We simply apply a first order implicit-explicit 

fn+1 _ fn T n+l 

1 +v-v x p = — (g[f n+i ] - r +i ), 



(IMEX) scheme, 
(2.2) 



At 

I f°(x,v) = f {x,v 



This can be written as 

p T- n+1 A/ 
(2-3) f n+1 = — -j— — [/" - Atv ■ V x f n ] -\ — =r^Q[f n+ \ 

y I J £ + T n+l At U J 1 £ + T n+l At U J ' 

where G{f n+l ) is the anisotropic Maxwellian distribution computed from f n+1 . Although 
(|2.3p appears nonlinearly implicit, since the computation of f n+1 requires the knowledge of 
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G[f n+1 ], it can be solved explicitly. Specifically, upon multiplying (j2.3|) by <p(v) denned by 

(j)(v) := (l,v,- 

and use the conservation properties of Q and the definition of Q[f] in (j 1 . 3 j) . (|1.4p . we define 
the macroscopic quantity U by U := (p,pu,E) computed from / and get [T2"| |2"5] 

r n+1 At 

<P(v)(r-AtvV x f n )dv + 



U 



n+l 



£ + T n +!At 



or simply 
(2.4) 



U 



n+l 



£ + T^At 



<f>(v) (f n -AtvV x f n )dv. 



cp(v)g[f n+1 ](v)dv, 



Thus U n+l can be obtained explicitly. This gives p n+l ,u n+l and T n+1 . Unfortunately, it is 
not enough to define C/[/ n+1 ] for which we need p n+l n+1 . Therefore, we define the tensor 
S by 



(2.5) 



-w+i 



v®vf n+1 dv = p n+1 (G n+1 + u n+1 ®u n+1 ) 
and multiply the scheme (j2.3|) by v ® v. Using the fact that 

v ® v Q[f] (v) dv = p (T + u g) it) , 



and (|2.5p . we get that 
(2.6) S n+1 = 



£ + (1 - i/)r n+1 At 



(1 _ j/Wn+l At 



At / v ®vv -V x f n dv 

n-1 ( T ™+ 1 H + u n+1 ®U n+l ) 



e+ (1 - v)T n+1 At 

Now Q[f n+1 ] can be obtained explicitly from U n+1 and S n+1 and then / n+1 from (12.3 



Finally the scheme reads 



(2.7) 



n+l 



u 



S n+1 



(u) (f n -AtvV x f n )dv, 



At 



£+ (1 - Z^T n+1 At 



(1 - I/) r n+1 At 
+ £+ (1 - v) T™+! At P 



v®vv ■ V x f n dv 
n+\ f T n+l u + u ™+l^ n «+l) 



/ 



n+l 



[r-At^-v.r] + 



r n+1 At 



Qlf 



n+ln 



e + T n+1 At L 1 £ + r n+1 At 

In summary, although (|2.2p is nonlinearly implicit, it can be solved explicitly, thus satisfies 
the second condition of an AP scheme. 

Now let us prove that the scheme (|2.7p preserves the correct asymptotic. 

Remark 2.1. When the IMEX scheme 12. 2\) is applied to the classical BGK equation il.6\) . 
U n+l in [2~4\ ) will completely defines M[P +1 ] given by < TOj) [Ql The steps after 

equation \2.J$ are new ideas introduced in this paper for the ES-BGK model. 



Proposition 2.2. Consider the numerical solution given by \2. 7|j. Then 
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(2.8) 



(i) For all e — > and At > 0, the distribution function f n+1 satisfies 

< f n+l (x,v) < max (||f|U \\G[f n+l ] |U) 
(ii) For all At > and f°, the distribution function f n converges to Ai[f n ], that is, 

\im f n = M[f n ] 

E— >0 

and the scheme gives a first order approximation in time of the compressible Euler 
system. 

(in) Moreover, if we asssume that \\f n — Ai n \\ = 0(e), for n > 2 and 

jjn+l _ jjn 



At 



< c, 



the scheme \2. 3\) asymptotically becomes a first order in time approximation of the 
compressible Navier-Stokes il.5\) given by 



p n+i _ p n 
At 



+ div x (p n u T 



0, 



p n + l u n+l - p n u n 
At 

E n+1 _ E r 



+ di Va; (p n u n ®u n + pTl) = ediv x (^(u"" 1 )) 



At 



+ div x ([E n + p n T n ] u n ) = ed\v x (p<j{u n - l )u n - 1 + kV x T n ~ l ) 



where r n p = p n / (1 — u) and r n n = (d v + 2)p n /2. 

Proof, (i) Let us first observe that / n+1 is a linear combination of f n (defined along charac- 
teristics) and G[f n ], thus we get the first assertion. 

To prove (ii), for any initial distribution function, we considet f n for n > 1 and compute 
the asymptotic limit of S n when e goes to zero in (|2.7p . it yields 



S n = p n pnj + n « ^ 



and using (|2.5|) . we also get 
thus 

rpn j 

Hence in the asymptotic limit e — > 0, the distribution function Q\f n \ becomes isotropic, which 
means that Q\f n \ converges to the classical Maxwellian 7W[/ n ]. 

Therefore, the solution at zeroth order is obtained by taking f n = M [f n ] in the conserva- 
tion laws (|2.4p . namely, 



U 



n+l 



4>{v) (M[f n ]-AtvV x M[f n ])dv, 



and the scheme for the macroscopic quantities reduces to the well-known approximation to 
the Euler equation 



(2.9) 



p n+l _ p n 
At 



+ div x (^n n ) = 0, 



p n + l u n+l - p n u n 
At 

E n+l _ E r, 



+ div x (p n u n ®u n + p n T n I) = 0, 



At 



+ div x ([E n + p n T n ] u n ) = 0. 
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Moreover, the temperature T n+1 satisfies the following 

J rpn+l _rnn J 

y At +fu n - V x T n + p n T n div x u n = 0(At). 

Now let us prove {in). The preservation of the asymptotic, that is the compressible Navier- 
Stokes equation, is based on the Chapman-Enskog method. It simply consists of expanding 
the distribution function f n into 



f n = M[f n ] + ef?, 



which implies that 



(2.10) 



fi(x,v)dv = 0, f /f (x, v) v dv = 0, / ff(x, v) \v | 2 dv = 0. 



These conditions are known as the compatibility conditions. Moreover, we also expand the 
stress tensor 



(2.11) 

and the heat flux 



9' 



(2.12) 
where 



Q n (x) 



v — u 



n\2 



(v -u n )f n {x,v)dv = + eQ?(x), 



©1=-/ f?(x,v)(y-u n )®(y-u n )dv, Q? 



1 7; Of 71 1 2 

f?(*,v) ] o (v-v»)dv 



and tr (0™ ) = 0. Inserting this latter expansion into the discrete conservation laws (j2.4j) . it 
gives 



p n+l _ p n 
At 



+ div x (p n u n ) = 0, 



p n+l u n+l _ p n u n 
At 

E n+l _ E n 



At 



+ div,. (p n u n ®u n + p n T n l) = -e div x (p n 9?) , 
+ div x {[E n + p n T n ] u n ) = -ediv x ( QJ + p n 9? u ) . 



For the application of the Chapman-Enskog method, the anisotropic Gaussian Q{ f) must be 
expanded with respect to e as well, 

Q\! n \ = M[f n ] + egl 

The next step is to insert these expansions into the scheme (|2.3f) for the ES-BGK eqution 
and use the compatibility conditions, and it yields for n > 1 

Mlf n ] — M\f n ~ 1 ~\ 
(2.13) U J ' [/ 1 + w • V^fr- 1 ] = r" (<tf - /?) 



At 



rn-l 



At J1 



One the one hand, the term is computed from G[f n ] as follows. First, the distribution 
function Q[f n ] contains the inverse matrix T n for which insertion of (|2.1ip yields 

r = ri + veQ\. 

Therefore, a Taylor expansion within terms of second order and using the fact that tr(0™) = 0, 
gives on the one hand 

det(T n ) = (T 1 *)*- + 0(e 2 ) 
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and on the other hand 



[TT 1 = I - + 0(e 



and therefore 



9i 



G[f n ]-M[f n ] 



j<r, 



(v-u n ye^ {v-u n ). 



e 2 (T n ) 2 

On the other hand, to obtain the first order expression for the distribution function we 
need to consider the terms of order zero in (|2.13p . which we write for convenience as 



/r = 9i 



M[f n ]-M[f 



n-l] 



+ —.V x M[f n - 1 ] + 0(e). 



The differential d.M[/] of the Maxwellian is given by 

dp f\v — u\ 2 d v \ dT u — v 



dM[f] = M[f]d\og(M[f]) = M[f] 



— + 

P 



2T 



2 T T 



and with the assumption (|2.8p . we obtain 
M[f n ]-M[f n - 1 } 



At 



M[f 



n-l] 



1 p n - p n ~ X 1 



n n-l 



At 



v — u 



n-l\2 



u n 1 — v u n — u n 1 



^n— 1 



At 



+ O(At) 



n rpn—1 



At 



and 



v-V x M[f n ~'}=M[f n - 1 ] 



v ■ V x p n 1 / | /' - a 



P 



,n-l 



n-l|2 



2T n_1 



2 



■n— 1 ZL n ~^ 



ipn—l 



J^n—l 



Gathering the two latter equalities, it yields up to the order At 
M[f n ] -M[f n - 1] 



+ 



At 

Mir- 1 ] 

ipn—l 

Mir- 1 ] 



-i + vV x M[f n - 1 ] = 
(u n ~ l - v) V,^ 1 -v) + (u"" 1 



\v — u 



n-l\2 



d v + 2 



2T n ~ l 



v x r 



n-l 



P 



n-l 



p n - p n ~ X 
At 



n—1 „,n— 1\ 



+ diy c (p n U 



Mif"- 1 ] {v - u"- 1 ) 



T 



+ Mir- 1 ] 



n-l 



,.n _ „n-l i 

.n-l V7 „,n— 1 i x V7 l„n—\rpn—l 



At 



v — u 



n-l\2 



2T 



n-l 



1 / T n T 



T 



•n—l 



n rpn—1 



At 



+ u"' 1 • V x T n ~' ) + — div.u 
d v 



.n-l 



The last three lines may vanish due to the conservation laws (|2.9p up to the order of e and 
At. Thus, the result for the first order contribution to the distribution function is 



f n 
Jl 



M[f 



n-l] 



+ V 



j-n rpn-1 

Mir 

2\T 



- v) V^u"" 1 (m"" 1 -v) + (u* 1 - 1 - v) 



\v — u 



n-l\2 



d v + 2 



2T 



n—l 



n-l 



n\2 



v - u n ) 6™ (v - u n ) + O(At) + 0{e). 



It is straightforward to show that fi satisifes the compatibility relations (|2,10p . Hence, we 
get for the stress tensor 0" 



n—l rpn—1 



r(v-u n )(g)(v-u n )dv 



p n-l jr? 



a(u n - L ) + vp n Q n x + O(At) + 0(e) 
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that is, 



n-1 
(1 — V)T n 



Then the heat flux Q™ is given by 

r |„ _ „n|2 H 4-9 n n ~ l 

Qi = / ' _ ' (^-n n )/f^ = -^L^^V.T"- 1 + O(Ai) + O(e). 
jRdu ^ It 

The constitutive relations are laws of Navier-Stokes and Fourier where 

1 p n - 1 , 4+2 p™- 1 
u = and k = 

^ 1 - f T n 2 T n 

are viscosity and thermal conductivity, respectively. Thus, a first order approximation with 
respect to At and e is given by 

p «ey = -/i^/- 1 ) and q? = -^r- 1 

The Prandtl number is related to the coefficient v of the ES-BGK model by 

d v + 2 /i 1 



Pr 



2 K 1-1/ 



□ 



3. Numerical simulations 

In this section, we give three numerical examples for the ES-BGK equation in different 
asymptotic regimes in order to check the performance (in stability and accuracy) of our 
methods. We have implemented the first order scheme (|2.7p for the approximation of the ES- 
BGK equation. A classical second order finite volume scheme with slope limiters is applied 
for the transport operator. We present two numerical tests for a ld x x 2d v model and finally 
a non stationary 2d x x 2d v model. 

We will compare the numerical solution to the ES-BGK equations with the one obtained 
for the full Boltzmann equation with Maxwellian molecules using a spectral approximation 
|17l [2] and the one obtained for the compressible Navier-Stokes system using a second order 
finite volume scheme. 

To this aim, we need to choose the right value for r such that the viscosity and heat 
conductivity computed from the asymptotic limit of the ES-BGK model is the same as the 
ones corresponding to the full Boltzmann equation. According to [10], the viscosity computed 
for the full Boltzmann equation is given by 

m ^ T 



3vr A 2 (5)' 

where ^(5) ~ 0.436. On the other hand, the viscosity computed from the ES-BGK model is 

1 p 
V = i > 

1 — V T 

where v = — 1. Thus, we choose r such that both viscosity are equal, which leads to 

1 V 37T , . , 7T 

t = ±7— = — = A 2 (5)p ~ 0.925 - p. 

2(i B (T) 2V2 ^ ' 2 H 
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fo(x, v) v dvdx = 0. 



3.1. Approximation of smooth solutions. For this numerical test, we consider the ES- 
BGK equation in dimension ld x x 2d v on the torus 

8f 1 

-+ + vVxf = -Q(f,f), xe[-i,i],vem. 2 
f(t = o) = f , 

with periodic boundary conditions in x. The operator Q(f) is given by Q(f) = r [G[f] — f] 
where the parameter r is chosen in order that the viscosity p matches perfectly with one 
obtained to the full Boltzmann operator for Maxwellian molecules, that is r = 0.9 tt p/2. 
Define p g and T g with respect to the initial data fo by 

Pg = \ I I fo(x,v)dvdx, and T g = — [ [ f (x,v) \v\ 2 dvdx 
1 J -i Jm 2 2 Pg J -i Jr 2 

and assume for simplicity that 

I/ 1 

Whenever f(t, x, v) is a smooth solution to the Boltzmann or ES-BGK equation with periodic 
boundary conditions, one has the global conservation laws for mass, momentum and energy. 
These conservation laws are then enough to uniquely determine the stationary state of the 
model : the normalized global Maxwellian distribution 

(3.14) M a (v) = exp (-^) , v eR 2 . 

Our goal here is to investigate numerically the long-time behavior of the solution / and 
to compare the solution with the asymptotic behavior of the solution to the compressible 
Navier-Stokes equations (CNS). If / is any reasonable solution of the ES-BGK equation, 
satisfying certain a priori bounds of compactness (in particular, ensuring that no kinetic 
energy is allowed to leak at large velocities), then it is expected that / does indeed converge 
to the global Maxwellian distribution Ai g as t goes to +oo. 

Recently, Desvillettes and Villani |13j . Guo and Strain [18] were interested in the study 
of rates of convergence for the full Boltzmann equation. Roughly speaking in [13j . the au- 
thors proved that if the solution to the Boltzmann equation is smooth enough then (with 
constructive bounds) 

\\f(t)-M g \\ = 0(t-°°), 

which means that the solution converges almost exponentially fast to the global equilibrium 
(namely with polynomial rate 0(t~ r ) with r as large as wanted). Moreover in [13], Desvillettes 
and Villani conjectured that time oscillations should occur on the evolution of the relative 
local entropy 



W-l(t) = j /log dxdv > 



where Mi is the local Maxwellian distribution in the sense that the constants p, u and T 
appearing there depend on time t and position x 

/„ .v/ \ P(t, x ) ( \v — u(t,x)\ 2 

(3.15) Mi(t,x,v) = / ' —s exp ' 



2irT(t,x) *\ 2T(t,x) 

In fact their proof does not rule out the possibility that the entropy production undergoes 
important oscillations in time, and actually most of the technical work is caused by this 
possibility. 

To estimate the speed of convergence to the global equilibrium and the possibility that 
oscillations also occur on the difference between global and local equilibria (since the global 
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relatve entropy is nonincreasing) , we prefer to investigate the behavior of the following quan- 
tity 

£(t) := \\p(t) - Pg \\ L i, 

which makes sense both for the solution to the Boltzmann or ES-BGK equation but also for 
the compressible Navier-Stokes equation. 

Here, we performed simulations to the full Boltzmann equation with a fast spectral method 
[IT] . to the ES-BGK equation with the scheme (|2.3p and to the compressible Navier-Stokes 
system with a WENO solver in a simplified geometry (one dimension in space, two dimension 
in velocity, periodic boundary conditions) to check numerically if such oscillations occur. 
Clearly this test is challenging for a numerical method due to the high accuracy required to 
capture such oscillating behavior at the kinetic regime. 

Then, we consider an initial datum as a perturbation of the global equilibrium M g 



fo(x,v) 



1 + Aq sin(7r x) 



exp 



\V - U \ 



+ exp 



\V + Uq\ 



,x £ [-1, l],v G 



2T ) l \ 2T 

with A§ = 0.5, To = 0.125 and no = (1/2, 1/2). We have chosen v = —1 such that the Prandtl 
number of the ES-BGK model corresponds to the Prandtl number of the 2d v Boltzmann 
operator, that is Pr = 0.5. 




0.001 



Time evolution of II p(t) - p 
5 



10 15 




Time evolution of II pit) -p g II ji 
5 10 



(c) 



(d) 



Figure 1. Influence of the Knudsen number e: distance between the local 
density p(t) and the global density at equilibrium p g using 100 x 32 x 32 for 
(a) e=0.5; (b) e = 0.1; (c) e = 0.05 and (d) e = 0.01. 



In Figure [H one can indeed observe oscillations on the quantitiy £{t) for the full Boltzmann 
equation, the ES-BGK model and also for the compressible Navier-Stokes system obtained 
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from the asymptotic of the ES-BGK equation (jl.9p . The strength of the oscillations does not 
depend on the Knudsen number e. The superimposed curves yield the time evolution of the 
£ (t) for t G [0,20]; the first plot corresponds to e = 0.5, the second one to e = 0.1, the third 
one e = 0.05 and the last one e = 0.01. 

On the one hand, for e = 0.5, which corresponds to a rarefied regime, the behavior of 
£ (t) strongly differs between the kinetic and hydrodynamic models. The results for ES-BGK 
and the full Boltzmann equations agree very well in this rarefied regime, which illustrates 
perfectly the efficiency and accuracy of the ES-BGK model. 

On the other hand, for smaller values of e, the different numerical approximations give 
roughly the same results and the ES-BGK model and the compressible Navier-Stokes system 
become very close. Finally for e ~ 0.01, the different kinetic (Boltzmann and ES-BGK) and 
hydrodynamic (Euler or Navier-Stokes) models agree very well. 

Further note that the equilibration is much more rapid when the Knudsen number e is 
large, and that the convergence seems to be exponential. 

3.2. The Riemann problem. This test deals with the numerical solution to the ld x x 
2d v ES-BGK equation. The operator Q(f) is given by Q(f) = t [£/[/] — /] where r = 
0.9 7rp/2. We present several numerical simulations for one dimensional Riemann problem, 
with different Knudsen numbers, from rarefied regime to the fluid regime. 

Here, the initial data corresponding to the ES-BGK equations are given by the isotropic 
Maxwellian distributions computed from the following macroscopic quantities 



with the Mach number M = 2.5. We perform several computations for e = 5 x 10 ,10 



We present a comparison between the numerical solution to the Boltzmann equation ob- 
tained using a spectral scheme [IT], the approximation to the compressible Navier-Stokes 
system obtained using a fifth order WENO and our first order implicit method (j2.3|) for the 
ES-BGK model with v = 1/2, for which the Prandtl number of the ES-BGK model is the 
same as the one for the 2d v Boltzmann operator. Let us note that the viscosity and conduc- 
tivity used for the numerical simumation of the compressible Navier-Stokes system are the 
ones obtained from the Chapman-Enskog expansion. 

In Figuer [21 we take e = 5 x 10 _1 and choose a time step At = 0.001 satisfying the CFL 
condition for the transport part (with n x = 200). For such a value of e, the problem is not 
stiff and this test is only performed to compare the accuracy of our scheme (j2.3j) with the 
different models. We present several snapshots of the density, mean velocity, temperature 
and heat flux 



at different time t = 0.1, 0.25 and 0.4. We observe that, for a short time t = 0.1, the 
numerical approximation of macroscopic quantities and heat flux given by (|2.3p for the ES- 
BGK model are relatively close to the numerical solution to the Boltzmann equation. The 
front speed and the shape of the temperature with two bumps are very well approximated 
by the ES-BGK model, and the heat flux given by the Boltzmann equation and the ES-BGK 
model are different from the ones given by the compressible Navier-Stokes system. Clearly, 
the compressible Navier-Stokes system is not adequate in this rarefied regime whereas the 
ES-BGK model gives very satisfying results. 

For a larger time t > 0.25, the distribution in velocity / to the ES-BGK model and to 
the Boltzmann equation agree very well and the macroscopic quantities are already well 
approximated by the solution to the compressible Navier-Stokes system. In Figure [3j we 




10" 3 . 
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represent the x — v x projection of f — M[f] at time t = 0.25 obtained for the Boltzmann 
equation and the ES-BGK model. The distribution function becomes particularly far from 
the equilibrium at the front of the shock in velocity and then propagates in the computational 
domain. As can be observed on the macroscopic quantities, the solution to the Boltmzann 
equation is very close to the one obtained from the ES-BGK model. 

Finally, at the kinetic regime our method (|2.3p gives the same accuracy as a standard first 
order fully explicit scheme for the ES-BGK model or full Boltzmann equation without any 
additional computational effort. Of course, the computational effort needed for the ES-BGK 
models is much smaller than the one for an accurate discretization of the full Boltzmann 
equation. 




(3) (4) 

Figure 2. Riemann problem (e = 5 x 10 _1 ): crosses (+) represent the nu- 
merical solution to the Boltzmann equation obtained with our method (j2.3|) . 
stars (x) represent the solution to the ES-BGK model and lines is the solution 
corresponding to the compressible Navier-Stokes system. Evolution of (1) the 
density p, (2) mean velocity u, (3) temperature T and (4) heat flux Q at time 
t = 0.1, 0.25 and 0.4. 

Now, we investigate the cases of small values of e for which an explicit scheme requires 
the time step to be of order 0{e). In order to evaluate the accuracy of our method (|2.3p in 
the Navier-Stokes regime (for small e <C 1 but not negligible), we compared the numerical 
solution for e = 10 _1 with one obtained by the approximation of the compressible Navier- 
Stokes system derived from the ES-BGK model since the viscosity and heat conductivity are 
in that case explicitly known [3]. 

Therefore, in Figure HI we report the numerical results for e = 10 _1 and make comparison 
between the numerical solution obtained with the scheme (|2.3p and the one obtained with a 
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~1i 

H °' 5 1 1 



0.15 
0.2 

0.25 



(1) 



(2) 



Figure 3. Riemann problem (e = 5 x 10 1 ): x — v x projection of / — A^[/] 
at time 0.25 for the (1) Boltzmann equation and (2) ES-BGK model. 



high order explicit method for the compressible Navier-Stokes. In this case, the behavior of 
macroscopic quantities (density, mean velocity, temperature and heat flux) agree very well 
even if the time step is at least ten times larger with our method (|2,3p . 

Finally in Figure [SJ we choose e = 1. 10~ 3 . In this problem, the density, mean velocity and 
temperature are relatively close to the one obtained with the approximation of the Navier- 
Stokes system. Even the qualitative behavior of the heat flux agrees well with the heat flux 
corresponding to the compressible Navier-Stokes system k,V x T, with k = pT (see Figure [5]), 
yet some differences can be observed, which means that the use of ES-BGK models to derive 
macroscopic models has a strong influence on the heat flux. 

3.3. Flow around a cylinder. This example has been considered in for instance [26J. The 
computational domain is set to be [—20, 20] x [—20, 20]. The cylinder is centered at the origin, 
with a diameter of 2. 

We consider an incoming flow at the boundary \\x\\ = 8 with the following conditions: 
Pi = 1 ; Ui = {My/2Xi^) T , Ti = 1 with M = 0.1 and M = 0.5. The freestream Knudsen 
number ranges from e = 0.1 to e = 1CP 3 . 

Concerning boundary condition, the wall of cylinder is considered as a pure diffusive 
boundary conditions 



f(t,x,v) 



p(t,x) 



cxp 



2T„ 



v ■ n x < 0, and 



1. 



where p is computed such that the global flux is zero at the boundary (and mass is preserved) 
and T w = 1.05. 

To start the calculation take a uniform initial solution equal to the values defined by the 
boundary conditions: 

|2^ 



fo(x,v) 



exp 



v G 



1 < llxll < 



(2ttT 1 )*'/ 2 ~" r V 2T, 
Then, we solve the kinetic equations for the different grid densities considered, until a steady 
state is reached. 

We define the Mach number from the macroscopic quantities, computing the moments of 
the distribution function with respect to v G M 2 , by 

IL.II2 



M 



7 T' 



where c := y^yT is the sound speed. 
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(3) (4) 

Figure 4. Riemann problem (e = 1CP 1 ): crosses (+) represent the numerical 
solution to the Boltzmann equation obtained with our method fj2.3[> . stars 
(x) represent the solution to the ES-BGK model and lines is the solution 
corresponding to the compressible Navier-Stokes system. Evolution of (1) the 
density p, (2) mean velocity u, (3) temperature T and (4) heat flux Q at time 
t = O.f, 0.25 and 0.4. 

We apply our numerical scheme (|2.3p to the ES-BGK equation and plot the numerical 
results in the following figures (Fig. I6|7l8p and the solution can be compared to the numerical 
solution of the Euler equations [26] . 

Figure [6] shows the density contours at different times for a free streaming Mach number 
M = 0.1 whereas the Knudsen number is e = 0.01. The reflecting shock, and the Mach shock 
can all be identified. However, the reflecting shock is not identified as kinetic, since both 
density and temperature are large on the shock, which makes its distribution much closer to 
the Maxwellian than the other two shocks. Also, the two seperation points, where the gas 
is the most rarefied, are well captured and correctly identified as kinetic. In Figure the 
contour plot of the local Mach number, the density and the temperature is shown when the 
steady state is reached. 

4. Conclusion 

In this paper we present an accurate deterministic method for the numerical approxi- 
mation of the space inhomogeneous, time dependent ES-BGK equation. The method is a 
temporally implicit-explicit scheme to deal with the stiffness of the collision operator. The 
computational cost of the implicit part is close to an explicit one, without using any nonlinear 
algebraic system solver, by utilizing the particular structure of the ES-BGK operator. This 
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Figure 5. Riemann problem (e = 1CP 3 ): crosses (+) represent the numerical 
solution to the Boltzmann equation obtained with our method fj2.3[> . stars 
(x) represent the solution to the ES-BGK model and lines is the solution 
corresponding to the compressible Navier-Stokes system. Evolution of (1) the 
density p, (2) mean velocity u, (3) temperature T and (4) heat flux Q at time 
t = O.f, 0.2 and 0.3. 



effective time discretization allows the treatment of problems with a broad range of mean 
free path. Moreover, the numerical results, and the comparison with other techniques, show 
the effectiveness of the present method for a wide class of problems. 
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